{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contact of an elastic membrane with a rigid indenter: an augmented-Lagrangian approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py:235: QuadratureRepresentationDeprecationWarning: \n",
      "*** ===================================================== ***\n",
      "*** FFC: quadrature representation is deprecated! It will ***\n",
      "*** likely be removed in 2018.1.0 release. Use uflacs     ***\n",
      "*** representation instead.                               ***\n",
      "*** ===================================================== ***\n",
      "  issue_deprecation_warning()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual: 1.04457776733e-05\n",
      "Residual: 8.09266832477e-06\n",
      "Residual: 6.445405364e-06\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-f238365fa626>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobstacle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m \u001b[0msolve_Uzawa\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-1-f238365fa626>\u001b[0m in \u001b[0;36msolve_Uzawa\u001b[0;34m(tol)\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     68\u001b[0m         \u001b[0mlocal_project\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mppos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m         \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mform\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mJ\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mJ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     70\u001b[0m         \u001b[0;31m#gap.assign(project(h-u, P))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     71\u001b[0m         \u001b[0;31m#p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() +\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    298\u001b[0m     \u001b[0;31m# tolerance)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    299\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclasses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEquation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 300\u001b[0;31m         \u001b[0m_solve_varproblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    301\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    302\u001b[0m     \u001b[0;31m# Default case, just call the wrapped C++ solve function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36m_solve_varproblem\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNonlinearVariationalSolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    348\u001b[0m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msolver_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 349\u001b[0;31m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    351\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "parameters[\"form_compiler\"][\"representation\"] = 'quadrature'\n",
    "import warnings\n",
    "from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning\n",
    "warnings.simplefilter(\"once\", QuadratureRepresentationDeprecationWarning)\n",
    "\n",
    "ppos = lambda x: (x-abs(x))/2.\n",
    "Heav = lambda x: (1-x/abs(x))/2.\n",
    "\n",
    "N = 40\n",
    "mesh = UnitSquareMesh(N, N, \"crossed\")\n",
    "\n",
    "obstacle = Expression(\"y+pow(x[0]-0.5,2)+pow(x[1]-0.5, 2)\", y=-0.001, degree=2)\n",
    "\n",
    "alpha = 100.\n",
    "Nitermax = 10\n",
    "\n",
    "V = FunctionSpace(mesh, \"CG\", 2)\n",
    "#P = FunctionSpace(mesh, \"Quadrature\", quadrature_scheme=\"default\")\n",
    "Pe = FiniteElement(\"Quadrature\", mesh.ufl_cell(), degree=2, quad_scheme='default')\n",
    "P = FunctionSpace(mesh, Pe)\n",
    "\n",
    "u = Function(V)\n",
    "u_old = Function(V)\n",
    "du = TrialFunction(V)\n",
    "u_ = TestFunction(V)\n",
    "h = Function(V)\n",
    "\n",
    "p = Function(P)\n",
    "p_old = Function(P)\n",
    "gap = Function(P)\n",
    "\n",
    "corr = Constant(0.)\n",
    "rho = Function(P)\n",
    "dxm = dx(metadata={\"quadrature_degree\": 2, \"quadrature_scheme\": \"default\"})\n",
    "lagrangian = 0.5*inner(grad(u), grad(u))*dx + p*ppos(h-u)*dx + alpha/2.*ppos(h-u)**2*dx\n",
    "form = inner(grad(u), grad(u_))*dx - p*Heav(h-u)*u_*dxm\n",
    "J = derivative(form, u, du)\n",
    "\n",
    "def local_project(v, V, u=None):\n",
    "    dv = TrialFunction(V)\n",
    "    v_ = TestFunction(V)\n",
    "    a_proj = inner(dv, v_)*dxm\n",
    "    b_proj = inner(v, v_)*dxm\n",
    "    solver = LocalSolver(a_proj, b_proj)\n",
    "    solver.factorize()\n",
    "    if u is None:\n",
    "        u = Function(V)\n",
    "        solver.solve_local_rhs(u)\n",
    "        return u\n",
    "    else:\n",
    "        solver.solve_local_rhs(u)\n",
    "        return\n",
    "\n",
    "def border(x, on_boundary):\n",
    "    return on_boundary\n",
    "bc = DirichletBC(V, Constant(0), border)\n",
    "\n",
    "def solve_Uzawa(tol=2e-6):\n",
    "    t = 1.\n",
    "    for i in range(Nitermax):\n",
    "        u_old.assign(u)\n",
    "        p_old.assign(p)\n",
    "        told = t\n",
    "\n",
    "        local_project(p + alpha*ppos(h-u), P, p)\n",
    "        solve(form == 0, u, bc, J=J)\n",
    "        #gap.assign(project(h-u, P))\n",
    "        #p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() + \n",
    "        #                           alpha*gap.vector().get_local())\n",
    "        #direction = assemble((h-u)*(p-p_old)*dx)\n",
    "        #t = (1+sqrt(1+4*told**2))/2.\n",
    "        #if direction >= 0:\n",
    "        #    corr=(t-1.)/told\n",
    "        #else:\n",
    "        #    print \"    Restart...\"\n",
    "        #    corr = 0\n",
    "        #    t = 1.\n",
    "        #rho.vector()[:] = (1+corr)*p.vector().get_local() - corr*p_old.vector().get_local()\n",
    "        res = max(0*errornorm(u, u_old), assemble((p- p_old)**2*dxm))\n",
    "        print \"Residual:\", res\n",
    "        if res < tol:\n",
    "            break\n",
    "\n",
    "h.interpolate(obstacle)     \n",
    "solve_Uzawa()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl = plot(u)\n",
    "plt.colorbar(pl)\n",
    "plt.show()\n",
    "\n",
    "pl = plot(local_project(p, V))\n",
    "plt.colorbar(pl)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "\n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-c6af918a3cf7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0mobstacle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobstacle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0msolve_Uzawa\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-12-d2d4361d090f>\u001b[0m in \u001b[0;36msolve_Uzawa\u001b[0;34m(tol)\u001b[0m\n\u001b[1;32m     43\u001b[0m         \u001b[0mtold\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m         \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mform\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m         \u001b[0;31m#gap.assign(project(h-u, P))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;31m#p.vector()[:] = np.minimum(0*p_old.vector().get_local(), rho.vector().get_local() +\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    298\u001b[0m     \u001b[0;31m# tolerance)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    299\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mufl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclasses\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mEquation\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 300\u001b[0;31m         \u001b[0m_solve_varproblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    301\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    302\u001b[0m     \u001b[0;31m# Default case, just call the wrapped C++ solve function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python2.7/dist-packages/dolfin/fem/solving.pyc\u001b[0m in \u001b[0;36m_solve_varproblem\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    347\u001b[0m         \u001b[0msolver\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNonlinearVariationalSolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    348\u001b[0m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msolver_parameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 349\u001b[0;31m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    350\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    351\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: \n\n*** -------------------------------------------------------------------------\n*** DOLFIN encountered an error. If you are not able to resolve this issue\n*** using the information listed below, you can ask for help at\n***\n***     fenics-support@googlegroups.com\n***\n*** Remember to include the error message listed below and, if possible,\n*** include a *minimal* running example to reproduce the error.\n***\n*** -------------------------------------------------------------------------\n*** Error:   Unable to solve nonlinear system with NewtonSolver.\n*** Reason:  Newton solver did not converge because maximum number of iterations reached.\n*** Where:   This error was encountered inside NewtonSolver.cpp.\n*** Process: 0\n*** \n*** DOLFIN version: 2017.2.0\n*** Git changeset:  unknown\n*** -------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "for d in np.linspace(1e-4, 1e-3, 11):\n",
    "    obstacle.y = d\n",
    "    h.interpolate(obstacle)\n",
    "    solve_Uzawa()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<dolfin.cpp.la.Vector; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Vector > *' at 0x7f401e7570c0> >"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assemble(form)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
